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ABSTRACT. We study systems of equations of the form X\ = fi(Xi, . . . , X n ), . . . , X„ = 
f n (X\, . . . , X n ) where each /; is a polynomial with nonnegative coefficients that add up to 1. The 
least nonnegative solution, say fj,, of such equation systems is central to problems from various areas, 
like physics, biology, computational linguistics and probabilistic program verification. We give a 
simple and strongly polynomial algorithm to decide whether fi = (1, . . . , 1) holds. Furthermore, we 
present an algorithm that computes reliable sequences of lower and upper bounds on converging 
linearly to fi. Our algorithm has these features despite using inexact arithmetic for efficiency. We 
report on experiments that show the performance of our algorithms. 



1. Introduction 



We study how to efficiently compute the least nonnegative solution of an equation system of 
the form 

Xi = fi (X\ , . . . , X n ) . . . X n = f n (Xi , . . . , X n ) , 
where, for every i G {1, . . . , n}, fi is a polynomial over X\, . . . , X n with positive rational coef- 
ficients that add up to The solutions are the fixed points of the function / : IR n — >• W 1 with 
/ = (/i, . . . , /„). We call / a probabilistic system of polynomials (short: PSP). E.g., the PSP 



f(Xi,X 2 ) 

induces the equation system 
Xx 



1 11 11 

-XiX 2 + - , ^2X2 + -Xx + - 



7^X1X2 + i X2 = \X2X2 + 4X1 + \ . 

Obviously, 1 = (1, . . . , 1) is a fixed point of every PSP. By Kleene's theorem, every PSP has a 
least nonnegative fixed point (called just least fixed point in what follows), given by the limit of the 
sequence 0, /(0), /(/(D)), . . . 

PSPs are important in different areas of the theory of stochastic processes and computational 
models. A fundamental result of the theory of branching processes, with numerous applications in 
physics, chemistry and biology (see e.g. [|9j |2l), states that extinction probabilities of species are 
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1 Later, we allow that the coefficients add up to at most 1. 
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equal to the least fixed point of a PSP. The same result has been recently shown for the probability 
of termination of certain probabilistic recursive programs (7J |6). The consistency of stochastic 
context-free grammars, a problem of interest in statistical natural language processing, also reduces 
to checking whether the least fixed point of a PSP equals 1 (see e.g. iTTTlD . 

Given a PSP / with least fixed point /U/, we study how to efficiently solve the following two 
problems: (1) decide whether [if = 1, and (2) given a rational number e > 0, compute lb, ub € Q n 
such that lb < [if < ub and ub — lb < I (where u < v for vectors u, v means < in all 
components). While the motivation for Problem (2) is clear (compute the probability of extinction 
with a given accuracy), the motivation for Problem (1) requires perhaps some explanation. In the 
case study of Section R31 we consider a family of PSPs, taken from [9], modelling the neutron 
branching process in a ball of radioactive material of radius D (the family is parameterized by D). 
The least fixed point is the probability that a neutron produced through spontaneous fission does 
not generate an infinite "progeny" through successive collisions with atoms of the ball; loosely 
speaking, this is the probability that the neutron does not generate a chain reaction and the ball does 
not explode. Since the number of atoms in the ball is very large, spontaneous fission produces many 
neutrons per second, and so even if the probability that a given neutron produces a chain reaction is 
very small, the ball will explode with large probability in a very short time. It is therefore important 
to determine the largest radius D at which the probability of no chain reaction is still 1 (usually 
called the critical radius). An algorithm for Problem (1) allows to compute the critical radius using 
binary search. A similar situation appears in the analysis of parameterized probabilistic programs. 
In (7J HI it is shown that the question whether a probabilistic program almost surely terminates can 
be reduced to Problem (1). Using binary search one can find the "critical" value of the parameter 
for which the program may not terminate any more. 

Etessami and Yannakakis show in Q that Problem (1) can be solved in polynomial time by 
a reduction to (exact) Linear Programming (LP), which is not known to be strongly polynomial. 
Our first result reduces Problem (1) to solving a system of linear equations, resulting in a strongly 
polynomial algorithm for Problem (1). The Maple library offers exact arithmetic solvers for LP and 
systems of linear equations, which we use to test the performance of our new algorithm. In the 
neutron branching process discussed above we obtain speed-ups of about one order of magnitude 
with respect to LP. 

The second result of the paper is, to the best of our knowledge, the first practical algorithm for 
Problem (2). Lower bounds for [if can be computed using Newton's method for approximating a 
root of the function f(X) — X. This has recently been investigated in detail Q [lOl El- However, 
Newton's method faces considerable numerical problems. Experiments show that naive use of exact 
arithmetic is inefficient, while floating-point computation leads to false results even for very small 
systems. For instance, the PReMo tool [12], which implements Newton's method with floating- 
point arithmetic for efficiency, reports [if > 1 for a PSP with only 7 variables and small coefficients, 
although [if < 1 is the case (see Section l3TI ). 

Our algorithm produces a sequence of guaranteed lower and upper bounds, both of which con- 
verge linearly to [if. Linear convergence means that, loosely speaking, the number of accurate bits 
of the bound is a linear function of the position of the bound in the sequence. The algorithm is 
based on the following idea. Newton's method is an iterative procedure that, given a current lower 
bound lb on [if, applies a certain operator Af to it, yielding anew, more precise lower bound Af(lh). 
Instead of computing A/"(lb) using exact arithmetic, our algorithm computes two consecutive New- 
ton steps, i.e., N(N(lh)), using inexact arithmetic. Then it checks if the result satisfies a carefully 
chosen condition. If so, the result is taken as the next lower bound. If not, then the precision is 
increased, and the computation redone. The condition is eventually satisfied, assuming the results 
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of computing with increased precision converge to the exact result. Usually, the repeated inexact 
computation is much faster than the exact one. At the same time, a careful (and rather delicate) 
analysis shows that the sequence of lower bounds converges linearly to 

Computing upper bounds is harder, and seemingly has not been considered in the literature be- 
fore. Similarly to the case of lower bounds, we apply / twice to ub, i.e., we compute /(/(ub)) with 
increasing precision until a condition holds. The sequence so obtained may not even converge to [it. 
So we need to introduce a further operation, after which we can then prove linear convergence. 

We test our algorithm on the neutron branching process. The time needed to obtain lower and 
upper bounds on the probability of no explosion with e = 0.0001 lies below the time needed to 
check, using exact LP, whether this probability is 1 or smaller than one. That is, in this case study 
our algorithm is faster, and provides more information. 

The rest of the paper is structured as follows. We give preliminary definitions and facts in 
Section [2] Sections [3] and |4]present our algorithms for solving Problems (1) and (2), and report on 
their performance on some case studies. Section [5] contains our conclusions. The full version of the 
paper, including all proofs, can be found in |4). 

2. Preliminaries 

Vectors and matrices. We use bold letters for designating (column) vectors, e.g. v G M n . We write 
s with s G M for the vector (s, . . . , s) T G W 1 (where T indicates transpose), if the dimension n 
is clear from the context. The z-th component of v G R n will be denoted by Vj. We write x = y 
(resp. x < y resp. x -< y) if x, = y» (resp. x« < y^ resp. x$ < y^) holds for alH G {1, . . . , n}. 
By x < y we mean x < y and x/y. By W nxn we denote the set of real matrices with m rows 
and n columns. We write Id for the identity matrix. For a square matrix A, we denote by p(A) 
the spectral radius of A, i.e., the maximum of the absolute values of the eigenvalues. A matrix is 
nonnegative if all its entries are nonnegative. A nonnegative matrix A G M. nxn is irreducible if for 
every k, I G {1, . . . , n} there exists an i G N so that (A 1 )^ ^ 0. 

Probabilistic Systems of Polynomials. We investigate equation systems of the form 

X\ = f\(Xi, . . . ,X n ) ... X n = f n (Xi, . . . , X n ), 

where the are polynomials in the variables X\, . . . , X n with positive real coefficients, and for 
every polynomial fa the sum of its coefficients is at most 1. The vector / := (/i, . . . , f n ) T is called 
a probabilistic system of polynomials (PSP for short) and is identified with its induced function 
/ : 1" -»• W 1 . If Xi, . . . ,X n are the formal variables of /, we define X := (X\, . . . , X n ) T 
and Var(/) := {X\, . . . , X n }. We assume that / is represented as a list of polynomials, and each 
polynomial is a list of its monomials. If S C {X\ , . . . , X n }, then fg denotes the result of removing 
the polynomial fi(Xi, . . . ,X n ) from / for every Xi S; further, given x G R n and B G M nxn , 
we denote by xg and Bss the vector and the matrix obtained from x and B by removing the entries 
with indices i such that Xi S. The coefficients are represented as fractions of positive integers. 
The size of / is the size of that representation. The degree of / is the maximum of the degrees of 
fi, . . . , f n . PSPs of degree (resp. 1 resp. >1) are called constant (resp. linear resp. superlinear). 
PSPs / where the degree of each /j is at least 2 are called purely superlinear. We write /' for the 
Jacobian of /, i.e., the matrix of first partial derivatives of /. 

Given a PSP /, a variable Xi depends directly on a variable Xj if Xj "occurs" in fi, more 
formally if J^- is not the constant 0. A variable Xi depends on Xj if Xi depends directly on Xj 
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or there is a variable such that X{ depends directly on and X^ depends on Xj. We often 
consider the strongly connected components (or SCCs for short) of the dependence relation. The 
SCCs of a PSP can be computed in linear time using e.g. Tarjan's algorithm. An SCC S of a PSP / 
is constant resp. linear resp. superlinear resp. purely superlinear if the PSP / has the respective 
property, where / is obtained by restricting / to the S-components and replacing all variables not 
in S by the constant 1. A PSP is an scPSP if it is not constant and consists of only one SCC. Notice 
that a PSP / is an scPSP if and only if /'(l) is irreducible. 

A fixed point of a PSP / is a vector x > with /(x) = x. By Kleene's theorem, there exists 
a least fixed point fit of f, i.e., [if < x holds for every fixed point x. Moreover, the sequence 
0) /(0)j /(/(0)), • • • converges to [if. Vectors x with x < /(x) (resp. x > /(x)) are called pre- 
fixed (resp. post-fixed) points. Notice that the vector 1 is always a post-fixed point of a PSP /, due 
to our assumption on the coefficients of a PSP. By Knaster-Tarski's theorem, [if is the least post- 
fixed point, so we always have < [if < 1. It is easy to detect and remove all components i with 
{[if)i = by a simple round-robin method (see e.g. Q), which needs linear time in the size of /. 
We therefore assume in the following that [if y 0. 

3. An algorithm for consistency of PSPs 

Recall that for applications like the neutron branching process it is crucial to know exactly 
whether [if = 1 holds. We say a PSP / is consistent if [if = 1; otherwise it is inconsistent. 
Similarly, we call a component i consistent if ([if)i = 1. We present a new algorithm for the 
consistency problem, i.e., the problem to check a PSP for consistency. 

It was proved in that consistency is checkable in polynomial time by reduction to Linear 
Programming (LP). We first observe that consistency of general PSPs can be reduced to consistency 
of scPSPs by computing the DAG of SCCs, and checking consistency SCC-wise Q: Take any 
bottom SCC S, and check the consistency of fs- (Notice that fs is either constant or an scPSP; 
if constant, fs is consistent iff fs = 1, if an scPSP, we can check its consistency by assumption.) 
If fs is inconsistent, then so is /, and we are done. If fs is consistent, then we remove every fi 
from / such that Xi G S, replace all variables of S in the remaining polynomials by the constant 1, 
and iterate (choose a new bottom SCC, etc.). Note that this algorithm processes each polynomial at 
most once, as every variable belongs to exactly one SCC. 

It remains to reduce the consistency problem for scPSPs to LP. The first step is: 

Proposition 3.1. 10171 An scPSP f is consistent iff p(f (1)) < 1 (i.e., iff the spectral radius of 
the Jacobi matrix f evaluated at the vector 1 is at most 1 ). 

The second step consists of observing that the matrix f'(\) of an scPSP / is irreducible and non- 
negative. It is shown in [7] that p(A) < 1 holds for an irreducible and nonnegative matrix A iff the 
system of inequalities 

Ax>x + T,x>0 (3.1) 
is infeasible. However, no strongly polynomial algorithm for LP is known, and we are not aware 
that (13.11 ) falls within any subclass solvable in strongly polynomial time (8|. 

We provide a very simple, strongly polynomial time algorithm to check whether p(f'(l)) < 1 
holds. We need some results from Perron-Frobenius theory (see e.g. Q). 

Lemma 3.2. Let A € R™*™ be nonnegative and irreducible. 

(1) p(A) is a simple eigenvalue of A. 

(2) There exists an eigenvector v >- with p(A) as eigenvalue. 
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(3) Every eigenvector v y has p(A) as eigenvalue. 

(4) For all a, f3 G R \ {0} and v > 0: ifav < Av < /3v, f/j<?n a < p(A) < /3. 
The following lemma is the key to the algorithm: 

Lemma 3.3. Let A G R nxn be nonnegative and irreducible. 

(a) Assume there is v£ff'\ {0} such that (Id - A)v = 0. Then p(A) < 1 iff v >- or v -< 0. 

(b) Assume v = is the only solution of (Id — A)v = 0. Then there exists a unique x G M n 
such that (Id — A)x = 1, and p(A) < 1 iff x > 1 awci Ax < x. 

P roo/ 

(a) From (Id — A)v = it follows Av = v. We see that v is an eigenvector of A with 
eigenvalue 1. So p(A) > 1. 

(<=): As both v and — v are eigenvectors of A with eigenvalue 1, we can assume w.l.o.g. 
that v y 0. By Lemma |3T2T 3), p(A) is the eigenvalue of v, and so p(A) = 1. 
(=>): Since p(A) < 1 and p(A) > 1, it follows that p(A) = 1. By Lemma 13111) and 
(2), the eigenspace of the eigenvalue 1 is one-dimensional and contains a vector x y 0. So 
v = a ■ x for some a G R, a / 0. If a > 0, we have v y 0, otherwise v -< 0. 

(b) With the assumption and basic facts from linear algebra it follows that (Id — A) has full 
rank and therefore (Id — A)x = 1 has a unique solution x. We still have to prove the second 
part of the conjunction: 

(<=): Follows directly from Lemma I3T21 4). 

(=^): Let p(A) < 1. Assume for a contradiction that p(A) = 1. Then, by Lemma I3T21T ). 
the matrix A would have an eigenvector v / with eigenvalue 1, so (Id — A)v = 0, 
contradicting the assumption. So we have, in fact, p(A) < 1. By standard matrix facts 
(see e.g. O), this implies that (Id — A) -1 = A* = X^o^* exists, and so we have 
x = (Id - A)- 1 ! = A*T > T. Furthermore, Ax = J2Zi A ^ < T,Zo A ^ = x - ■ 

In order to check whether p(A) < 1, we first solve the system (Id — A)v = using Gaussian 
elimination. If we find a vector v ^ such that (Id— A)\ = 0, we apply Lemma[33Ja). If v = is 
the only solution of (Id — A)v = 0, we solve (Id — A)v = 1 using Gaussian elimination again, and 
apply Lemma l3"3l b). Since Gaussian elimination of a rational n-dimensional linear equation system 
can be carried out in strongly polynomial time using 0(n 3 ) arithmetic operations (see e.g. O), we 
obtain: 

Proposition 3.4. Given a nonnegative irreducible matrix A G M. nxn , one can decide in strongly 
polynomial time, using 0(n 3 ) arithmetic operations, whether p(A) < 1. 

Combining Propositions 13.11 and 13.41 directly yields an algorithm for checking the consistency 
of scPSPs. Extending it to multiple SCCs as above, we get: 

Theorem 3.5. Let f(Xi, . . . ,X n ) be a PSP. There is a strongly polynomial time algorithm that 
uses 0(n 3 ) arithmetic operations and determines the consistency of f. 

3.1. Case study: A family of "almost consistent" PSPs 

In this section, we illustrate some issues faced by algorithms that solve the consistency problem. 
Consider the following family of scPSPs, n > 2: 

= ( 0.5Xf + 0.1X* + 0.4 , COIA^ 2 + 0.5X 2 + 0.49 , . . . , 0.01X 2 _ 1 + 0.5X n + 0.49 ) T . 
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n = 25 


n = 100 


n = 200 


n = 400 


n = 600 


n = 1000 


Exact LP 


< 1 sec 


2 sec 


8 sec 


67 sec 


208 sec 


> 2h 


Our algorithm 


< 1 sec 


< 1 sec 


1 sec 


4 sec 


10 sec 


29 sec 



Table 1: Consistency checks for /i( n )-systems: Runtimes of different approaches. 



It is not hard to show that /i (n) (p) ~< p holds for p = (1 - 0.02 n , . . . , 1 - 0.02 2n ~ 1 ) T , so we have 
A^M ^ 1 by Proposition 14.41 i.e., the are inconsistent. 

The tool PReMo [12] relies on Java's floating-point arithmetic to compute approximations of 
the least fixed point of a PSP We invoked PReMo for computing approximants of /i fo ( n ) for different 
values of n between 5 and 100. Due to its fixed precision, PReMo's approximations for ^ h ( n ) are 
> 1 in all components if n > 7. This might lead to the wrong conclusion that is consistent. 

Recall that the consistency problem can be solved by checking the feasibility of the system (13.11 ) 
with A = /'(l). We checked it with lp_solve, a well-known LP tool using hardware floating-point 
arithmetic. The tool wrongly states that (13.11) has no solution for -systems with n > 10. This is 
due to the fact that the solutions cannot be represented adequately using machine number precision^ 
Finally, we also checked feasibility with Maple's Simplex package, which uses exact arithmetic, and 
compared its performance with the implementation, also in Maple, of our consistency algorithm. Ta- 
ble [Qshows the results. Our algorithm clearly outperforms the LP approach. For more experiments 
see Section 1431 

4. Approximating fi / with inexact arithmetic 

It is shown in [7 ] that /i j may not be representable by roots, so one can only approximate fij. In 
this section we present an algorithm that computes two sequences, (lb^)j and (ub^)j, such that 
lbW < Hf < ubW and Hindoo ub^ - lb^ = 0. In words: lb^ and ub^ are lower and upper 
bounds on [if, respectively, and the sequences converge to/if. Moreover, they converge linearly, 
meaning that the number of accurate bits of lb^ and ub^* 1 are linear functions of i. (The number of 
accurate bits of a vector x is defined as the greatest number k such that \ {fJ-f — x)j |/| (nf)j [ < 2~ k 
holds for all j € {1, . . . ,n}.) These properties are guaranteed even though our algorithm uses 
inexact arithmetic: Our algorithm detects numerical problems due to rounding errors, recovers from 
them, and increases the precision of the arithmetic as needed. Increasing the precision dynamically 
is, e.g., supported by the GMP library HI. 

Let us make precise what we mean by increasing the precision. Consider an elementary op- 
eration g, like multiplication, subtraction, etc., that operates on two input numbers x and y. We 
can compute g(x,y) with increasing precision if there is a procedure that on input x, y outputs a 
sequence (x,y), (x,y), . . . that converges to g(x,y). Note that there are no requirements 
on the convergence speed of this procedure — in particular, we do not require that there is an i 
with g^(x,y) = g(x,y). This procedure, which we assume exists, allows to implement floating 
assignments of the form 

z *~ g(x, y) such that <f>(z) 
with the following semantics: z is assigned the value g( l \x, y), where i > 1 is the smallest index 
such that <p(g^ (x, y)) holds. We say that the assignment is valid if 4>(g(x, y)) holds and <f> involves 

2 The mentioned problems of PReMo and lp.solve are not due to the fact that the coefficients of cannot be 
properly represented using basis 2: The problems persist if one replaces the coefficients of h^"' by similar numbers 
exactly representable by machine numbers. 
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only continuous functions and strict inequalities. Our assumption on the arithmetic guarantees that 
(the computation underlying) a valid floating assignment terminates. As "syntactic sugar", more 
complex operations (e.g., linear equation solving) are also allowed in floating assignments, because 
they can be decomposed into elementary operations. 

We feel that any implementation of arbitrary precision arithmetic should satisfy our require- 
ment that the computed values converge to the exact result. For instance, the documentation of 
the GMP library 12 states: "Each function is defined to calculate with 'infinite precision' followed 
by a truncation to the destination precision, but of course the work done is only what's needed to 
determine a result under that definition." 

To approximate the least fixed point of a PSP, we first transform it into a certain normal form. A 
purely superlinear PSP / is called perfectly superlinear if every variable depends directly on itself 
and every superlinear SCC is purely superlinear. The following proposition states that any PSP / 
can be made perfectly superlinear. 

Proposition 4.1. Let f be a PSP of size s. We can compute in time 0(n ■ s) a perfectly superlinear 
PSP f with Var(f) = Var(f) U {X} of size 0(n ■ s) such that p,f = (nf)var(f)- 

4.1. The algorithm 

The algorithm receives as input a perfectly superlinear PSP / and an error bound e > 0, and 
returns vectors lb,ub such that lb < fit < ub and ub — lb < e. A first initialization step 
requires to compute a vector x with -< x -< /(x), i.e., a "strict" pre-frxed point. This is done in 
Section l4.1.1l The algorithm itself is described in Section 14.1.21 

4.1.1. Computing a strict pre-fixed point. Algorithm [TJcomputes a strict pre-fixed point: 

Algorithm 1: Procedure computeStrictPref ix 
Input: perfectly superlinear PSP / 
Output: x with -< x -< /(x) -< T 
x 4- 0; 

while / xdo 

Z <- {i | 1 < i < n,/i(x) = 0}; 
P<- 01 1 < i < n,/i(x) > 0}; 
yz <- 0; 

yp *~ /p(x) such that -< y P -< f P (y) -< 1; 
_ x i- y; 



Proposition 4.2. Algorithm\l\is correct and terminates after at most n iterations. 

The reader may wonder why Algorithm Q] uses a floating assignment yp ^~ /p(x), given that 
it must also perform exact comparisons to obtain the sets Z and P and to decide exactly whether 
yp ~< /p(y) holds in the such that clause of the floating assignment. The reason is that, while we 
perform such operations exactly, we do not want to use the result of exact computations as input for 
other computations, as this easily leads to an explosion in the required precision. For instance, the 
size of the exact result of fp(y) may be larger than the size of y, while an approximation of smaller 
size may already satisfy the such that clause. In order to emphasize this, we never store the result 
of an exact numerical computation in a variable. 
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4.1.2. Computing lower and upper bounds. Algorifhm[T]uses Kleene iteration 0, /(0), /(/(0)), . . . 
to compute a strict pre-fixed point. One could, in principle, use the same scheme to compute lower 
bounds of /ij, as this sequence converges to [if from below by Kleene's theorem. However, conver- 
gence of Kleene iteration is generally slow. It is shown in [7] that for the 1-dimensional PSP / with 
f(X) = 0.5X 2 + 0.5 we have [if = 1, and the i-th Kleene approximant kW satisfies /«W < 1 — i. 
Hence, Kleene iteration may converge only logarithmically, i.e., the number of accurate bits is a 
logarithmic function of the number of iterations. 

In it was suggested to use Newton's method for faster convergence. In order to see how 
Newton's method can be used, observe that instead of computing [if, one can equivalently compute 
the least nonnegative zero of f(X) — X. Given an approximant x of [if, Newton's method first 
computes g^(X), the first-order linearization of / at the point x: 

< 7 W(X) = /(x) + /'(x)(X-x) 

The next Newton approximant y is obtained by solving X = g^(X), i.e., 

y = X+(/d-/'(x))- 1 (/(x)-x) . 

We write A/V(x) := x + (Id — /'(x)) _1 (/(x) — x), and usually drop the subscript of J\[f. If 
< [if is any pre-fixed point of /, for instance = 0, we can define a Newton sequence 
(i) )i by setting = A/"0 W ) for i > 0. It has been shown in HSSIIH that Newton sequences 

converge at least linearly to [if. Moreover, we have < i/W < f(y®) < [if for all i. 

These facts were shown only for Newton sequences that are computed exactly, i.e., without 
rounding errors. Unfortunately, Newton approximants are hard to compute exactly: Since each 
iteration requires to solve a linear equation system whose coefficients depend on the results of the 
previous iteration, the size of the Newton approximants easily explodes. Therefore, we wish to 
use inexact arithmetic, but without losing the good properties of Newton's method (reliable lower 
bounds, linear convergence). 

Algorithm [2] accomplishes these goals, and additionally computes post-fixed points ub of /, 
which are upper bounds on [if. Let us describe the algorithm in some detail. The lower 
bounds are stored in the variable lb. The first value of lb is not simply 0, but is computed by 
computeStrictPref ix(/), in order to guarantee the validity of the following floating assign- 
ments. We use Newton's method for improving the lower bounds because it converges fast (at least 
linearly) when performed exactly. In each iteration of the algorithm, two Newton steps are per- 
formed using inexact arithmetic. The intention is that two inexact Newton steps should improve the 
lower bound at least as much as one exact Newton step. While this may sound like a vague hope 
for small rounding errors, it can be rigorously proved thanks to the such that clause of the floating 
assignment in line [4] The proof involves two steps. The first step is to prove that M(M(\b)) is a 
(strict) post-fixed point of the function g(X) = /(lb) + f'(lb)(X - lb), i.e., N(Af(lb)) satisfies 
the first inequality in the such that clause. For the second step, recall that A/"(lb) is the least fixed 
point of g. By Knaster-Tarski's theorem, .A/" (lb) is actually the least post-fixed point of g. So, our 
value x, the inexact version of J\f(J\f (lb)), satisfies x > .A/" (lb), and hence two inexact Newton 
steps are in fact at least as "fast" as one exact Newton step. Thus, the lb converge linearly to [if. 

The upper bounds ub are post-fixed points, i.e., /(ub) < ub is an invariant of the algorithm. 
The algorithm computes the sets Z and P so that inexact arithmetic is only applied to the compo- 
nents i with /t(ub) < 1. In the P-components, the function / is applied to ub in order to improve 
the upper bound. In fact, / is applied twice in line |9l analogously to applying M twice in line 0] 
Here, the such that clause makes sure that the progress towards [if is at least as fast as the progress 
of one exact application of / would be. One can show that this leads to linear convergence to [if. 
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Algorithm 2: Procedure calcBounds 



Input: perfectly superlinear PSP /, error bound e > 

Output: vectors lb, ub such that lb < fj,t < ub and ub — lb < e 

1 lb 4— computeStrictPref ix(/); 

2 ub T; 

3 while ub — lb ^ e do 
x ^M(M{lb)) such that /(lb) + /'(lb)(x - lb) ■< x ■< /(x) -< I; 
lb^x; 

Z <-{*|l<f<n,/i(ub) = l}; 
P 4- {i | 1 < i < n,/i(ub) < 1}; 

yz <-T; 

yp *~ /p(/(ub)) such that / P (y) -< y P -< f P (ub); 
forall superlinear SCCs S of / with y$ = 1 do 
t <- T - lb 5 ; 
if f' ss (T)t >- t then 



ub 



y; 



minjgg(/^ g (l)t - t)j 
2 • max ie 5(/ s (2)) i 



t such that f s (y) -< ys H 1; 



The rest of the algorithm (lines [T0lfT3l) deals with the problem that, given a post-fixed ub, the 
sequence ub, /(ub), /(/(ub)), . . . does not necessarily converge to \xj. For instance, if f(X) = 
0.75X 2 + 0.25, then fi f = 1/3, but 1 = /(l) = /(/(l)) = Therefore, the if-statement 
of Algorithm [2] allows to improve the upper bound from 1 to a post-fixed point less than 1, by 
exploiting the lower bounds lb. This is illustrated in Figure Q] for a 2-dimensional scPSP /. The 
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(a) (b) 
Figure 1: Computation of a post-fixed point less than 1. 
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dotted lines indicate the curve of the points (X\,X2) satisfying X\ = 0. 8X1X2 + 0.2 and X2 = 
OAXl +O.IX2 + 0.5. Notice that fif ~< T = /(T). In Figured] (a) the shaded area consists of 
those points lb where /'(1)(1 — lb) >- 1 — lb holds, i.e., the condition of line [12] One can show 
that \Xf must lie in the shaded area, so by continuity, any sequence converging to fit, in particular 
the sequence of lower bounds lb, finally reaches the shaded area. In Figure [T] (a) this is indicated 
by the points with the square shape. Figure Q] (b) shows how to exploit such a point lb to compute 
a post-fixed point ub -< 1 (post-fixed points are shaded in Figure [T](b)): The post-fixed point ub 
(diamond shape) is obtained by starting at 1 and moving a little bit along the straight line between 
1 and lb, cf. line [13] The sequence ub, /(ub), /(/(ub)), . . . now converges linearly to jif. 

Theorem 4.3. Algorithm\2\terminates and computes vectors lb, ub such that lb < jij < ub and 

ub — lb < e. Moreover, the sequences of lower and upper bounds computed by the algorithm both 
converge linearly to fj, f. 

Notice that Theorem l4.3l is about the convergence speed of the approximants, not about the time 
needed to compute them. To analyse the computation time, one would need stronger requirements 
on how floating assignments are performed. 

The lower and upper bounds computed by Algorithm [2] have a special feature: they satisfy 
lb -< /(lb) and ub > /(ub). The following proposition guarantees that such points are in fact 
lower and upper bounds. 

Proposition 4.4. Let f be a perfectly superlinear PSP. Let < x < 1. If~x.< /( x )> then x -< nj. 
Ifx > /(x), then x > p,t. 

So a user of Algorithm[2]can immediately verify that the computed bounds are correct. To summa- 
rize, Algorithm [2] computes provably and even verifiably correct lower and upper bounds, although 
exact computation is restricted to detecting numerical problems. See Section FOI for experiments. 

4.2. Proving consistency using the inexact algorithm 

In Section [3] we presented a simple and efficient algorithm to check the consistency of a PSP. 
Algorithm [2] is aimed at approximating /ij, but note that it can also prove the inconsistency of a 
PSP: when the algorithm sets ubj < 1, we know < 1. This raises the question whether 

Algorithm [2] can also be used for proving consistency. The answer is yes, and the procedure is 
based on the following proposition. 

Proposition 4.5. Let f be an scPSP. Let t >- be a vector with /'(l)t < t. Then f is consistent. 

Proposition 14 . 5 1 can be used to identify consistent components. 

Use Algorithm [2] with some (small) e to compute ub and lb. Take any bottom SCC S. 

• If /'(1)(1 — lbs) < 1 — lbs, mark all variables in S as consistent and remove the S- 
components from /. In the remaining components, replace all variables in S with 1. 

• Otherwise, remove S and all other variables that depend on S from /. 
Repeat with the new bottom SCC until all SCCs are processed. 

There is no guarantee that this method detects all i with (fif)i = 1. 
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n 


20 
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50 


100 


20 
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50 


100 


20 
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50 


100 


20 


10 

50 


100 
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Cons, check (Alg. Sec© 


< 1 


< 1 


2 


< l 


< l 


2 


< l 


< l 


2 


< i 


< l 


2 


Cons, check (exact LP) 


< 1 


20 


258 


< l 


22 


124 


< l 


16 


168 


< i 


37 


222 


Approx. Q D (e = 10~ 3 ) 


< 1 


< 1 


4 


2 


8 


32 


l 


5 


21 


l 


4 


17 


Approx. Q D {e = 10" 4 ) 


< 1 


< 1 


4 


2 


8 


34 


2 


7 


28 


i 


6 


23 



Table 2: Runtime in seconds of various algorithms on different values of D and n. 



4.3. Case study: A neutron branching process 

One of the main applications of the theory of branching processes is the modelling of cascade 
creation of particles in physics. We study a problem described by Harris in |9]. Consider a ball of 
fissionable radioactive material of radius D. Spontaneous fission of an atom can liberate a neutron, 
whose collision with another atom can produce further neutrons etc. If D is very small, most 
neutrons leave the ball without colliding. If D is very large, then nearly all neutrons eventually 
collide, and the probability that the neutron's progeny never dies is large. A well-known result shows 
that, loosely speaking, the population of a process that does not go extinct grows exponentially over 
time with large probability. Therefore, the neutron's progeny never dying out actually means that 
after a (very) short time all the material is fissioned, which amounts to a nuclear explosion. The 
task is to compute the largest value of D for which the probability of extinction of a neutron born 
at the centre of the ball is still 1 (if the probability is 1 at the centre, then it is 1 everywhere). This 
is often called the critical radius. Notice that, since the number of atoms that undergo spontaneous 
fission is large (some hundreds per second for the critical radius of plutonium), if the probability of 
extinction lies only slightly below 1, there is already a large probability of a chain reaction. Assume 
that a neutron born at distance £ from the centre leaves the ball without colliding with probability 
Z(£), and collides with an atom at distance r) from the centre with probability density 77). Let 
further f(x) = ^2 i> QPi.x t , where pi is the probability that a collision generates i neutrons. For a 
neutron's progeny to go extinct, the neutron must either leave the ball without colliding, or collide 
at some distance 77 from the centre, but in such a way that the progeny of all generated neutrons goes 
extinct. So the extinction probability Qd{Q of a neutron born at distance £ from the centre is given 
by g), p. 86: 

Qd(0 = K0+ I R&v)HQd(v)) dv 

Jo 

Harris takes f(x) = 0.025 + 0.830x + 0.07x 2 + 0.05x 3 + 0.025x 4 , and gives expressions for both 
and rj). By discretizing the interval [0, D] into n segments and replacing the integral by 
a finite sum we obtain a PSP of dimension n + 1 over the variables {Qn(jD /n) [ < j < n}. 
Notice that Qd(0) is the probability that a neutron born in the centre does not cause an explosion. 

Results. For our experiments we used three different discretizations n = 20, 50, 100. We applied 
our consistency algorithm from Section|3]and Maple's Simplex to check inconsistency, i.e., to check 
whether an explosion occurs. The results are given in the first 3 rows of Table|2j Again our algorithm 
dominates the LP approach, although the polynomials are much denser than in the -systems. 

We also implemented Algorithm [2] using Maple for computing lower and upper bounds 
on Qd(0) with two different values of the error bound e. The runtime is given in the last two 
rows. By setting the Digits variable in Maple we controlled the precision of Maple's software 
floating-point numbers for the floating assignments. In all cases starting with the standard value 
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of 10, Algorithm [2] increased Digits at most twice by 5, resulting in a maximal Digits value of 20. 
We mention that Algorithm [2] computed an upper bound -< 1, and thus proved inconsistency, after 
the first few iterations in all investigated cases, almost as fast as the algorithm from Section[3] 

Computing approximations for the critical radius. After computing Q rj (0) for various values of D 
one can suspect that the critical radius, i.e., the smallest value of D for which (Jd(0) = 1, lies 
somewhere between 2.7 and 3. We combined binary search with the consistency algorithm from 
Section |3]to determine the critical radius up to an error of 0.01. During the binary search, the algo- 
rithm from Section|3]has to analyze PSPs that come closer and closer to the verge of (in)consistency. 
For the last (and most expensive) binary search step that decreases the interval to 0.01, our algorithm 
took <1, 1, 3, 8 seconds for n = 20, 50, 100, 150, respectively. For n = 150, we found the critical 
radius to be in the interval [2.981, 2.991]. Harris [9] estimates 2.9. 

5. Conclusions 

We have presented a new, simple, and efficient algorithm for checking the consistency of PSPs, 
which outperforms the previously existing LP-based method. We have also described the first al- 
gorithm that computes reliable lower and upper bounds on /xj. The sequence of bounds converges 
linearly to fif. To achieve these properties without sacrificing efficiency, we use a novel combina- 
tion of exact and inexact (floating-point) arithmetic. Experiments on PSPs from concrete branching 
processes confirm the practicality of our approach. The results raise the question whether our com- 
bination of exact and inexact arithmetic could be transferred to other computational problems. 

Acknowledgments. We thank several anonymous referees for pointing out inaccuracies and helping 
us clarify certain aspects of the paper. The second author was supported by the DFG Graduiertenkol- 
leg 1480 (PUMA). We also thank Andreas Reuss for proofreading the manuscript. 

References 

[1] GMP library, http://gmplib.org. 

[2] K. B. Athreya and P. E. Ney. Branching Processes. Springer, 1972. 

[3] A. Berman and R. J. Plemmons. Nonnegative Matrices in the Mathematical Sciences. SIAM, 1994. 

[4] J. Esparza, A. Gaiser, and S. Kiefer. Computing least fixed points of probabilistic systems of polynomials. Technical 
report, Technische Universitat Miinchen, Institut fur Informatik, 2009. 

[5] J. Esparza, S. Kiefer, and M. Luttenberger. Convergence thresholds of Newton's method for monotone polynomial 
equations. In Proceedings ofSTACS, pages 289-300, 2008. 

[6] J. Esparza, A. Kucera, and R. Mayr. Model checking probabilistic pushdown automata. In LICS 2004, pages 12-21. 
IEEE Computer Society, 2004. 

[7] K. Etessami and M. Yannakakis. Recursive Markov chains, stochastic grammars, and monotone systems of nonlin- 
ear equations. Journal of the ACM, 56(1): 1-66, 2009. 

[8] M. Grotschel, L. Lovasz, and A. Schrijver. Geometric Algorithms and Combinatorial Optimization. Springer, 1993. 

[9] T. E. Harris. The theory of branching processes. Springer, Berlin, 1963. 
[10] S. Kiefer, M. Luttenberger, and J. Esparza. On the convergence of Newton's method for monotone systems of 

polynomial equations. In Proceedings ofSTOC, pages 217-226. ACM, 2007. 
[11] C. D. Manning and H. Schuetze. Foundations of Statistical Natural Language Processing. MIT Press, June 1999. 
[12] D. Wojtczak and K. Etessami. PReMo: an analyzer for probabilistic recursive models. In TACAS, volume 4424 of 
Lecture Notes in Computer Science, pages 66-71. Springer, 2007. 



This work is licensed under the Creative Commons Attribution-NoDerivs License. To view a 
copy of this license, visit jhttp : //creativecommons . org/licenses/by-nd/3 . 0/| 



